#include "sets.h"
#include "utils.h"
#include <cmath>
using namespace std;

typedef pair<double,double> pdbl;

double norma2(const pdbl& p1, const pdbl& p2) {
	return sqr(p1.first-p2.first)+sqr(p1.second-p2.second);
}

vector<pair<bset,bomba> > getSetsT(const vector< pair<pol,pol> >& ptsP, double t, double r, double R) {
	int n=ptsP.size();
	vector< pdbl > pts(n);
	vector<pair<bset,bomba> > res;
	double r2=sqr(r);
	forn(i,n) pts[i]=pdbl(eval(ptsP[i].first, t), eval(ptsP[i].second, t));
	forn(j,n)forn(i,j+1) { //tomo las parejas en ambos sentidos para tener los 2 circulos posibles
		pdbl p=pdbl(pts[i].first-pts[j].first,pts[i].second-pts[j].second);
		double d2=sqr(p.first)+sqr(p.second);
		if (0.25*d2>r2) continue; //demasiado lejos
		double alpha=atan2(p.second,p.first); if (alpha<0) alpha+=2.0*PI;
		double x=sqrt(d2)/2.0;
		forn(s,2) {
			double y=sqrt(r2-0.25*d2)*(1-s*2); //probar raiz positiva y negativa
			pdbl q(x*cos(alpha)-y*sin(alpha),x*sin(alpha)+y*cos(alpha));
			q.first+=pts[j].first;
			q.second+=pts[j].second;
			//el circulo tiene centro q y radio r
			if (sqr(q.first)+sqr(q.second)<=sqr(R)) continue; //demasiado cerca!
			bset st;
			forn(k,n) if (norma2(q,pts[k])<=r2+eps2) st[k]=true;
			res.push_back(make_pair(st,(bomba){q.first,q.second,t}));
		}
	}
	return res;
}

vector<double> _cc;
bool comp(const pair<bset,bomba>& p1, const pair<bset,bomba>& p2) {
	double r=1.0/log(p1.second.t)-1.0/log(p2.second.t);
	forn(i,_cc.size()) if (p1.first[i] != p2.first[i]) 
		r+=p1.first[i]?_cc[i]:-_cc[i];
	return r>0.0;
}

bool compt(const pair<bset,bomba>& p1, const pair<bset,bomba>& p2) {
	return p1.second.t>p2.second.t;
}

void fix(vector<pair<bset,bomba> >& v, int csets, int m, double alpha) {
	//eliminar sets incluidos en otros
	sort(v.begin(),v.end(),compt);
	forn(i,v.size()) {
		forn(j,v.size())if (i!=j) {
			if ((v[j].first | v[i].first) == v[j].first) { //si v[i] incluido en v[j]
				swap(v[i],v[v.size()-1]);
				v.pop_back();
				i--;
				break;
			}
		}
	}
	if (csets>=v.size()) return;
	//rankear y quedarse con los csets mejores
	_cc=vector<double>(m);
	forn(i,v.size())forn(j,m) if (v[i].first[j]) _cc[j]++;
	forn(j,m) _cc[j]=alpha+1.0/_cc[j];
	sort(v.begin(), v.end(), comp);
	v.resize(csets);
}

vector<pair<bset,bomba> > getSets(const vector< pair<pol,pol> >& ptsP, 
	double r, double it, double et, int disc, int iter, int csets,
	double alpha, double R ) {
		r-=eps*2.0;
	
	double delta=(et-it)/(disc-1);
	
	vector<double> ts(disc),nts(disc);
	vector<pair<double,int> > ptsTs(disc);
	ts[0]=it; 
	forn(i,disc-1) ts[i+1]=ts[i]+delta;

	vector<pair<bset,bomba> > res;

	forn(itNr,iter) {
		forn(i,disc) {
			vector<pair<bset,bomba> > rr=getSetsT(ptsP,ts[i],r,R);
			ptsTs[i].first=0;
			ptsTs[i].second=i;
			forn(j,rr.size()) ptsTs[i].first+=sqr(rr[j].first.count()); 
			res.resize(res.size()+rr.size());
			forn(j,rr.size()) res[res.size()-rr.size()+j]=rr[j];
			if (res.size()>csets) fix(res,csets,ptsP.size(),alpha);
		}
		fix(res,csets,ptsP.size(),alpha);
		sort(ptsTs.begin(),ptsTs.end(), greater<pair<double,int> >());
		delta/=3.0;
		forn(i,disc) {
			nts[i]=ts[ptsTs[i].second];
			if (i&1) nts[i]-=delta; else nts[i]+=delta;
			if (nts[i]<it+delta) nts[i]=it+delta;
			if (nts[i]>et-delta) nts[i]=et-delta;
		}
		swap(nts,ts);
	}
	
	return res;
}
